Conductance of correlated systems: real-time dynamics in finite systems 
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Numerical time evolution of transport states using time dependent Density Matrix Renormaliza- 
tion Group (td-DMRG) methods has turned out to be a powerful tool to calculate the linear and 
finite bias conductance of interacting impurity systems coupled to non-interacting one-dimensional 
leads. Several models, including the Interacting Resonant Level Model (IRLM), the Single Impurity 
Anderson Model (SIAM), as well as models with different multi site structures, have been subject 
of investigations in this context. In this work we give an overview of the different numerical ap- 
proaches that have been successfully applied to the problem and go into considerable detail when we 
comment on the techniques that have been used to obtain the full I-V-characteristics for the IRLM. 
Using a model of spinless fermions consisting of an extended interacting nanostructure attached to 
non-interacting leads, we explain the method we use to obtain the current-voltage characteristics 
and discuss the finite size effects that have to be taken into account. We report results for the linear 
and finite bias conductance through a seven site structure with weak and strong nearest-neighbor 
interactions. Comparison with exact diagonalisation results in the non-interacting limit serve as a 
verification of the accuracy of our approach. Finally we discuss the possibility of effectively enlarging 
the finite system by applying damped boundaries and give an estimate of the effective system size 
and accuracy that can be expected in this case. 

PACS numbers: 73.63.-b, 72.10.Bg, 71.27.+a, 73.63.Kv 



OVERVIEW 



During the past decade improved experimental tech- 
niques have made the production of and measurements 
on one-dimensional systems possible [l| , and hence led to 
an increased theoretical interest in these systems. How- 
ever, the description of non-equilibrium transport prop- 
erties, like the finite bias conductance of an interacting 
nanostructure attached to leads, is a challenging task. 
In general, for non-interacting particles, the conductance 
can be extracted from the transmission of the single par- 
ticle levels 0-I3l ■ For interacting particles in small or low- 
dimensional structures where the screening of electrons is 
reduced, electron-electron correlations can no longer be 
neglected. Recently several methods to calculate the zero 
bias conductance of strongly interacting nanostructures 
have been developed. One class of approaches consists 
in extracting the conductance from an easier to calcu- 
late equilibrium quantity, e.g. the conductance can be ex- 
tracted from a persistent current calculation from 
phase shifts in NRG calculations [l^l, or from approx- 
imations based on the tunneling density of states [11 



Alternatively one can evaluate the Kubo formula within 
Monte-Carlo simulations [l2|, or from DMRG calcula- 
tions [l3l - fl5| . Linear conductance has also been inves- 
tigated using Functional Renormalization Group stud- 
ies (T6j . or by diagonalizing small clusters and attaching 
them to leads via a Dyson equation [T7| . 



In contrast, there are only a few methods available 
to get rigorous results for the finite bias conductance. 
While the problem has been formally solved by Meir 
and Wingreen using Keldysh Greens functions [I^, the 
evaluation of these formulas for interacting systems is 
generally based on approximations such as real time 
Keldysh RG [l^. Within the framework of time de- 
pendent density functional theory (td-DFT) and Keldysh 
Greens functions Stefanucci and Almbladh [l^, [2l| dis- 
cuss the extraction of conduction from real time simu- 
lations. The restriction to finite sized systems for cal- 
culating transport within td-DFT was also discussed by 
Di Vcntra and Todorov In [2^ Bushong, Sai, and 

Di Vcntra discuss the extraction of a finite bias current 
similar as discussed below in the framework of td-DFT. 
Weiss, Eckel, Thorwart and Egger Q discuss an itera- 
tive method based on the summation of real-time path 
integrals (ISPI) in order to address quantum transport 
problems out of equilibrium. Han and Heary [25j discuss 
strongly correlated transport in the Kondo regime using 
imaginary time Quantum Monte Carlo techniques. 

In this work we review the concept of calculating the 
finite bias conductance of nanostructures based on real 
time simulations [26l - l4l| within the framework of the 
DMRG j43 - |46| . It provides a unified description of 
strong and weak interactions and works in the linear and 
finite bias regime, as long as finite size effects are treated 
properly. The method was successfully applied to obtain 
results for the finite bias conductance in the interact- 
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FIG. 2: Exponential damping in the leads with tk = A~''^'^t. 
In the damping region, the hopping parameter is reduced by 
powers of the damping constant < K~^^^ < 1, while it is 
at the constant value t where connected to the nanostructure 
and at the constant value A~"''^t on the boundaries. 



FIG. 1: Interacting nanostructure • attached to non- 
interacting leads O (finite interaction Uc with the first lead 
site • allowed) and schematic density profile (green solid line) 
of the A^-particle wavepacket at initial time T = 0. The den- 
sity profile corresponds to the A'^-particle ground state of the 
Hamiltonian H + Hsd, cf. Eq. (O, where the bieis voltage 
enters as a local chemical potential Vsd (black dotted line). 



ing resonant level model, showing perfect agreement with 
analytical methods based on the Bethe ansatz [sl] . I-V- 
characteristics have been obtained for the single- impurity 
Anderson model using the adaptive td-DMRG-method 
[S^ l . Finite size effects and especially the impact of the 
possible combinations of tight binding leads with an even 
or odd number of sites coupled to the structure have been 
studied in detail in [sl] for a single impurity and for three 
quantum dots. Here, we show that finite size effects can 
be directly related to the structure of the single particle 
energy levels in non-interacting systems. 

In a first approach of time dependent dynamics within 
DMRG, Cazalilla and Marston integrated the time- 
dependent Schrodinger equation in the Hilbert space ob- 
tained in a finite lattice ground state DMRG calcula- 
tion [2^ . Since this approach docs not include the density 
matrix for the time evolved states, its applicability is very 
limited. Luo, Xiang and Wang [23| improved the method 
by extending the density matrix with the contributions of 
the wave function at intermediate time steps, restricting 
themselves to the infinite lattice algorithm. Schmitteck- 
ert [lOl showed that the calculations can be considerably 
improved by replacing the integration of the time depen- 
dent Schrodinger equation with the evaluation of the time 
evolution operator using a Krylov subspace method for 
matrix exponentials and by using the full finite lattice 
algorithm. 

An alternative approach is based on wave function pre- 
diction [131. There one first calculates an initial state 
with a static DMRG. One iteratively evolves this state 
by combining the wave function prediction with a time 
evolution scheme. In contrast to the above mentioned 
full td-DMRG, one only keeps the wave functions for two 
time steps in each DMRG step. Different time evolution 
schemes have been implemented in the past using ap- 



proximations like the Trotter decomposition [28l . |29| . |32| , 
or the Runge-Kutta method [3l| . Schneider and Schmit- 
teckert (40l.l48j combined the idea of the adaptive DMRG 
method with direct evaluation of the time evolution op- 
erator via a matrix exponential using Krylov techniques 
as described in Ref. [30|. Therefore the method involves 
no Trotter approximations, the time evolution is unitary 
by construction, and it can be applied to models beyond 
nearest-neighbor hopping. 

Concerning finite size effects, damped boundary con- 
ditions have been applied in order to obtain an increased 
effective system size in the regime of small bias voltage 
[sgI Isrj , where an improved scheme for linear conduc- 
tance was presented in |14{ . In the non-interacting case 
this can be traced back to a shift of the discrete single 
particle energy levels of the system towards the center 
of the cosine band. We demonstrate that this procedure 
can also be used when applying bias voltage of the order 
of magnitude of the band width when handled carefully. 



II. THE SYSTEM 

The Hamiltonian for the nanostructure is given by (S: 
the structure itself, L: leads, C: contacts) 
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where fij = c^jCj. Individual sites are labeled according to 
Fig. [1] A/pot = m — n is the size of the interacting nanos- 
tructure, Vg denotes a local external potential, which can 
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FIG. 3: Different initial conditions, corresponding to (a) 
Hinit. = H + FsD(iVL - iVR)/2 and (b) Hi„it. = H. The 
band width for the cosine band is it. Assuming a single par- 
ticle picture, we understand that in case (a), increasing the 
bias voltage Vsd to a value greater than the band width qual- 
itatively does not change the initial state, since all particles 
populate only one of the two leads, while for case (b), quench- 
ing the leads to different energies at the initial time prevents 
some particles (holes) from tunneling from one lead to the 
other because of energy conservation. For this reason there is 
no current flow in the extreme case of Vsd > 4t, cf. Fig. [S] 



be applied to the nanostructure, [/ is a nearest-neighbor 
interaction inside the nanostructure, and Uc is a nearest- 
neighbor interaction with the first lead sites. The hop- 
ping elements in the leads, the structure, and coupling of 
the structure to the leads are tj , tg , and tc , respectively. 
The hopping parameter in the leads tj is not necessarily 
constant to allow for the inclusion of damped boundary- 
conditions. This can be used to divide the leads in three 
areas. Fig. [5J here, two regions with constant hopping 
matrix element t and A~"/^t are smoothly coupled via a 
region of exponential damped hopping, which allows for 
increasing the resolution of the level spacing of the single 
particle energy levels on the energy scale A~"/^t. For 
hard-wall boundary-conditions, however, tj=t = const. 

The current operator Ij at an arbitrary bond j can 
be derived from the charge operator Qj = —ehj using a 
continuity equation. For the tight-binding Hamiltonian 
([T]) the current operator and its expectation value take 
the form 
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We define the current through the nanostructure as an 
average over the current in the left and right contacts to 
the nanostructure 



I{T) = [/„_i(T) + (r)]/2. 
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FIG. 4: Time dependent current through a single impurity 
coupled to noninteracting ID leads for vanishing gate volt- 
age Vg = 0. The system consists of M lattice sites and TV 
particles at nominal filling N/M = 0.5. We find three time 
domains: 1. an initial transient regime with decaying oscilla- 
tions, 2. a pseudo stationary current plateau and 3. finite size 
reflections, (a) Shortly after the initial switching of the bias 
voltage the time dependent behavior is dominated by oscilla- 
tions which decay to a constant current plateau on the time 
scale Ts (here: tc = 0.3t, M = 120). (b) The flnite size of the 
system leads to reflections at the boundaries. A wave packet 
that runs through the system starting at the impurity will be 
reflected at the boundaries and returns to the impurity after 
time Tr. This results in the typical pattern with recurrent 
sign changes of the current (here: tc = 0-5t, M — 60). 



III. INITIAL CONDITIONS AND TIME 
EVOLUTION 

Following the prescription implemented in (sol . l39j we 
add an external bias potential, namely the charge oper- 
ator, 



Hi 



SD 




(7) 



to the unperturbed Hamiltonian H and take the ground 
state l^-o) = |*(T = 0)) of H + Hsd, obtained by 
a standard finite lattice DMRG calculation, as initial 
state at time T = d^]. The minimization of the 
energy of the system leads to a charge imbalance in 
the right (source) and the left (drain) lead correspond- 
ing to Vsd, as sketched in Fig. [Sfa). Alternatively, the 
bias voltage also can be added to the time evolution. 
The initial state |\l/o) then has to be obtained as the 
ground state of the unperturbed Hamiltonian H, while 
the time evolution is performed using H + Hgu , cf. also 
Fig. [DJb). Starting from |^'o), the time evolution of the 
system results from the time evolution operator U{T) 
with |*(r)) = [/(T)|*o), which leads to flow of the ex- 
tended wave packet through the whole system until it 
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FIG. 5: I-V-characteristics for the resonant level model with 
tc = OAt and Uc = 0. The linear conductance is 1. The plot 
shows results for two different time evolution schemes: (a) 
the initial state I'I'o) of the system is the ground state of the 
Hamiltonian i/ + Vsd(-^Vl — iVR)/2, while the time evolution is 
performed as |*(T)) = exp(-i_&T)|'I'o). (b) the initial state 
l^'o) of the system is the ground state of the Hamiltonian H, 
while the time evolution is performed as |»I'(T)) — exp[— i(_ff + 
Vsn{N-L - 7V'R)/2)r]|*o). For values of the bias voltage much 
smaller than the band width the both approaches agree nicely. 
However, we find strong deviations when band edge effects 
come into play. Note that (a) corresponds to the situation 
of wide band metallic leads. Since our emphasis lies on the 
description of nanostructures attached to metallic leads we 
prefer to work in this approach. When describing situations 
with band gap materials as leads one should refer to approach 
(b). For further discussion see Fig. [3] and the text. 
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FIG. 6: Time dependent current through a single impurity 
coupled to noninteracting ID leads with tc = OAt and Uc = 
2.0t for different values of Vsd and vanishing gate voltage 
Vg = 0. The system consists of M = 48 lattice sites and 
TV particles at nominal filling N/M — 0.5. The current is 
obtained from a td-DMRG calculation by performing the time 
evolution on an initial non equilibrium state, using a DMRG 
projection scheme with a variable number of kept states 100 < 
A'^cut < 5600 with the discarded entropy kept below a 
certain value (here: Sd ^ 10~^; cf. also Fig. [7]). (a) The initial 
state l^o) corresponds to the situation sketched in Fig. |3ja) 
where |*I'o) is obtained as the ground state of //init. = H + 
Vsn{N-L - TV'r)/2, (c) The initial state |*o) is obtained as 
the ground state of Hinit L „ „• The current plateau 
we are looking for can be obtained more reliable when using 
prescription (a). 



is reflected at the hard wall boundaries as described in 
[30| . Corresponding to the two different schemes intro- 
duced before, C/ is given as either (a) tj{T) = Q-^HT/h 
(b) U{T) = c-HH+ffsD)T/;i^ 

The sudden switching of the bias voltage results in a 
ringing of the current in a transient time regime [49j , see 
also Fig. IHa). Here we show the short time behavior 
of the current through a single impurity coupled to two 
leads in a system with M = 120 lattice sites in total. This 
transient behavior with its characteristic oscillations de- 
cays on the time scale Tg oc F, where F is the width of 
the conductance peak. By smearing out the voltage drop 
over a few lattice one may reduce the influence of large 
momentum states. Furthermore, the finite size of the 
system leads to reflection of wave packets at the bound- 
aries, cf. Fig. |4l^b). A wave packet travelling with Fermi 
velocity vp from the impurity towards the boundaries 
will return to the impurity after a transit time given by 
Tr oc M/vf, which is the characteristic time scale for 
finite size effects appearing in the expectation value of 
time dependent observables. 

To compare the approaches (a) and (b), we show cur- 
rent voltage-characteristics in Fig.[5]for the resonant level 
model with a single impurity (Afoot = m — n = 1, ci. 
Fig. [T]) coupled to two leads via the hopping matrix el- 



ement tc = 0.4t and the gate voltage as well as the in- 
teraction set to Uc = Vg = 0. The dots correspond 
to results obtained numerically using exact diagonalisa- 
tion, while the lines correspond to analytic calculations 
included for comparison. Here, the straight line shows 
the current assuming linear scaling with Vsd with lin- 
ear conductance 5 = 1, while the curved line overlaid by 
the numerical results for approach (a) has been obtained 
using the Landauer-Biittiker approach, taking cosine- 
dispersion into account. 

The procedure of extracting the current from the nu- 
merical data will be described in the next section. Here 
we want to emphasize the different results we get for 
the I-V-curve for the two different cases. For the tight 
binding Hamiltonian the dispersion relation is given by 
ek = — 2icosfc, with a finite band width At. For the 
approach (a) this leads in the non-interacting case to a 
saturation of /(Vsd) for all values of the bias voltage 
Vsd > 4t. Further increasing Vsd beyond the band edge 
docs not change the initial occupation of energy levels. In 
contrast, for the case (b), the particles will be distributed 
equally over the left and the right lead in the initial state 
4*0 ), whereas the voltage enters in the time evolution op- 
erator. For small values of Vsd we find a good agreement 
for /(Vsd) for (a) and (b), while for Vsd }t 2i there is 
a mismatch which finds its expression in a current max- 
imum for < Vsd < 4t with a subsequent break down 
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FIG. 7: Maximum dimension A^cut of the DMRG projection 
scheme for an I-V-calculation necessary to keep the discarded 
entropy Sd below a certain value (here: Sd ^ 10~^) for dif- 
ferent configurations I to IV and different values of the bias 
voltage VsD, where we used 100 < A'cut < 5600 states as a sec- 
ond limitation. Here, the current through the contact links 
to a single impurity with tc = OAt is obtained for 70 time 
steps (AT = OAh/t) in a system with M = 48 lattices sites 
at half filling, (a) The initial state |^'o) is the ground state 
of -ffinit. = H + V'sd(A^l - Nb)/2, (c) I^o) is obtained as the 
ground state of Hinit. I, „ „. 



figure caption of Fig. [6l For both approaches (a) and 
(c), we find a time regime of (quasi) constant current. 
However, approach (a) has several advantages over (c): 
the current plateau is more consistent, which simplifies 
analysis, and to keep the discarded entropy Sd in the 
td-DMRG calculation below a predefined threshold, the 
number of states, which have to be kept in the DMRG, is 
considerably higher for (c) when compared to (a), mak- 
ing approach (c) computationally much more expensive. 
The latter point is illustrated in Fig. [71 where we compare 
the maximum dimension N^^t of the DMRG projection 
scheme that is necessary to keep Sd ^ 10"'^, for different 
values of the bias voltage Vsd, of the gate voltage Vg and 
of the interaction Uc- We always find a much smaller 
value of TVcut for (a) as compared to (c). 

Another problem of approach (c) is the discretization 
of the I-V-curve into steps resulting from the discrete 
single particle energy levels of the initial state. This could 
probably be handled using a procedure similar to the one 
described in section IVBl 

For these reasons we will use approach (a) throughout 
the remainder of this paper. 



IV. DIFFERENTIAL AND LINEAR 
CONDUCTANCE 



to / = for Vsd > This behavior has been predicted 
in [soj and can be understood from Fig. [3]^b), which ex- 
plains how energy conservation prevents particles (holes) 
to tunnel from one lead to the other which removes con- 
tributions to the current. (sEj. More recently, a detailed 
analysis of the negative differential conductance for the 
situation (b) has been carried out [5l|. In this work, it 
has been realised that the density of states in the leads 
adds a major contribution to the breakdown of the cur- 
rent. 

Moreover, there are other approaches to how the ini- 
tial state and the time evolution can be defined. For 
example, in addition to prescription (a), the coupling tc 
and the interaction Uc can be set to zero for the calcu- 
lation of l^'o). In this case (c), both leads as well as the 
structure are totally independent systems, and there is a 
very intuitive connection of Vsd and the difference of the 
particle number in the left and the right lead, because 
the isolated leads can be described in a single particle 
picture. The drawback of this approach, which adds a 
sudden switching of tc and Uc in addition to the switch- 
ing of Vsd at initial time T = 0, is an enhanced transient 
regime and therefore a reduced plateau of constant cur- 
rent that we need to extract the I-V-curve from. In Fig. [6] 
we compare the time dependent current obtained using 
the different initial conditions (a) and (c) for a single 
impurity coupled to two leads via tc = OAt, including 
a finite density-density interaction Uc = 2M, for differ- 
ent values of Vsd. To evaluate the time evolution of a 
system with finite interaction numerically, we used the 
td-DMRG method, with parameters as described in the 



For the calculation of the DC-conductance through the 
nanostructure the time evolution has to be carried out 
for sufficiently long times until a quasi-stationary state is 
reached and the steady state current / can be calculated. 
If the stationary state corresponds to a well-defined ap- 
plied external potential Vsd , the differential conductance 
is given by ^(Vsd) = edI{VsT))/dVsT)- In the limit of a 
small applied potential, Vsd ~^ Oj the linear conductance 
is given by (/(Vsd) = e/(VsD)/VsD. 

To discuss the general behavior of the time evolution 
from an initial nonequilibrium state we first consider the 
most simple case we can think of: transport through a 
single impurity. The current rises from zero and set- 
tles into a quasi-stationary state, Fig. [DJa). After the 
wavepackets have traveled to the boundaries of the sys- 
tem and back to the nanostructure, the current falls back 
to zero and changes sign, cf. Fig.[4jb). Additionally there 
is a third type of finite size oscillations. Fig. [8l Here we 
show the time dependent current for different configura- 
tions, from the leads to the impurity on a single (left or 
right) contact link, and through the impurity as defined 
in Eq. After the initial oscillations have decayed on 
the time scale Ts, the current through a single contact 
link shows remaining oscillations, with an amplitude de- 
pending on Vsd and Vg, and proportional to the inverse 
of the system size 1/M. The latter is demonstrated in 
Fig. [31 The period of the oscillation depends on the ap- 
plied bias voltage [compare Fig. [8l (b, c)] but is inde- 
pendent of the system size [Fig. [51 (b-d)] and of the gate 
potential [Fig. [TU], and is given by Tj = 27r?i/|VsD|. In 
the resonant tunneling case [Fig. [S^a), V^ = 0], the os- 
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FIG. 8: Current through a single impurity with tc = 0.3t 
at nominal filling N/M — 0.5 obtained from exact numeri- 
cal diagonalization (a-c), or DMRG including interaction (d), 
respectively - (a) for different system sizes M at bias volt- 
age VsD ~ 0-U and gate voltage Vg — 0. The black dashed 
line corresponds to the mean value of the fit values I for the 
left and right contact link, for AI = 60 lattice sites. The fit 
interval has to be chosen carefully - initial oscillations from 
the bias voltage switching and the finite transit time have to 
be taken into account. Even though the period of the finite 
size oscillations considerably exceed the system size M = 60 
for Vsu ~ 0.1, the fit current I is in nice agreement with the 
current plateau of the M = 120 system. However, finite size 
effects still have to be addressed (b, Vg = 0.3t, Vsd = O.lt, 
and c, = 0.3t, Vsd = 0.4t) since in general the fit cur- 
rent can strongly depend on the system size - in particular, 
a non-zero gate voltage changes the particle number density 
in the leads when the overall particle number is fixed. The 
same fit procedure can be applied to interacting systems (d, 
Uc = 2.0t, Vsd = OAt, Vg = 0.3t). 



dilations on the left and the right contact link cancel in 
the current average Eq. ^ due to a different sign in the 
amplitude of the oscillations /j, which does not hold in 
general [Fig. [51[b-d), 7^ 0], where the amplitude of the 
oscillations as a function of the gate potential Vg varies 
differently on the individual contact links, Fig. [TUl 
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FIG. 9: Oscillation amplitude 7j from fits as shown in Fig.[8l 
as function of the inverse system size 1/M for different values 
of Vsd , of the time dependent current through a single contact 
link to a single impurity, with tc ~ 0.5t and = 0. 



In Fig. \TU\ we plot the fit of the oscillation frequency 
cjj = 27r/Tj as a function of the gate potential Vg for a 
fixed value of Vsd, where we find ujj to be independent 
of the gate potential. To be precise, the fit nicely con- 
firms the above relation of Vsd and oscillation period. 
This periodic contribution to the current is reminiscent 
of the Josephson contribution in the tunneling Hamilto- 
nian, obtained by gauge transforming the voltage into a 
time dependent coupling ic{T) = ic e'^^°^/'* [52|. Like 
in a tunnel barrier in a superconductor, we have a phase 
coherent quantum system, namely the ground state at 
zero temperature. Instead of the superconducting gap 
we have a finite size gap resulting from the finite nature 
of the leads. Therefore the amplitude of this residual wig- 
gling vanishes proportional to the finite size gap provided 
by the leads. 

The stationary current is given by a fit to / + 
/j cos{2ttT/Tj + (p) with the fit-parameters tagged by a 
tilde, because the oscillation period Tj is known. In gen- 
eral, the density in the leads, and therefore also the cur- 
rent, depends on the system size and a finite size analysis 
has to be carried out in order to extract quantitative re- 
sults [Fig.[8](b,c), see also discussion of Fig. [18]. Only in 
special cases (symmetry, half filled leads, and zero gate 
potential) is the stationary current independent of the 
system size [Fig. [8] (a)]. 



V. FINITE SIZE EFFECTS 

Finite size effects such as the finite transit time of a 
wave packet traveling through the system and the peri- 
odic contribution to the current make it difficult to ob- 
tain a pseudo-stationary state where a constant current 
can be extracted from the time evolution. This prob- 
lem can be treated by a fit procedure as discussed in 
the previous section. However, in the small bias regime, 
where the amplitude of the oscillations is bigger than the 
(expected) current and the oscillation time Tj exceeds 
the transit time, this approach is unreliable. In section 
IVIII we discuss the possibility of effectively enlarging the 
system using damped boundary conditions (DBG) while 
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FIG. 10: Fit of the oscillation frequency ojj = 27r/r,j of the 
Josephson oscillations in a system with M = 120 lattice sites 
and a single resonant level with tc = 0.3t at a bias voltage 
VsD = QAt. The oscillation period extracted from the time 
evolution of the current is in excellent agreement with the 
analytical expression ojj = |VsD|/ft (dashed black line). The 
kinks that appear in ojj can be traced back to the fact that 
the amplitude of the oscillations Jj vanishes for Vg « iVsD/S 
at either the left or the right contact link. Then a fit of iDj 
does not work. The residual wiggling (its amplitude as well as 
its frequency) depends on the size and the position of the fit 
interval [Tminjlmax], and is therefore consistent with a finite 
fitting interval in time domain. Enlarging the fit intervall 
in conjunction with the system size reduces this effect (not 
shown here). 



keeping the system size M constant (cf. Fig. [2]). Further- 
more, the time evolution of the current strongly depends 
on the number of lattice sites of the leads being even 
or odd, Figs. [Til HSl In Fig- El wc compare this effect 
for a non-interacting two-dot structure for different sys- 
tem sizes in the regime of very small voltage Vsd ^ ^, 
where we consider three qualitatively different cases, (a) 
Tr < Tj, (b) Tr « Tj and (c) Tr > Tj, where Tr,Tj dc- 
note the transit time and oscillation period respectively, 
as discussed in Sec lIIII Since the number of single par- 
ticle energy levels is equal to the number of lattice sites, 
these relations are connected to Vsd and the level spacing 
Ae as, (a) Ae > ^sd, (b) Ae « Fsd and (c) Ae < 1/sd- 
Intuitively one would expect that the level discretisation 
must be small compared to the energy scales of interest, 
and indeed we find, that on the time scale T < Tr the nu- 
merical simulation fits best with the analytic result /lb 
obtained from the Landauer-Biittiker approach in case 
(c) (see FigfTTj). However, in all cases, the time evolution 
of the current depends on the different configurations of 
the leads with even or odd number of lattice sites. Two 
aspects must be distinguished: (1) the qualitative dif- 
ference in the time evolution depending on wether the 
number of lead sites is equal (as for the e2e and the o2o 
configuration), or unequal (as for the e2o and the o2e 
configuration), is clearly demonstrated in the figure. For 
the two-dot structure, this holds true even for Tr ^ Tj, 
Fig. [TT] (c). For the o2o and the e2e configurations we 



find a behavior where the current suddenly increases by a 
factor of ~ 2 after the transit time Tr, as opposed to the 
"expected" behavior with a sign change, seen for the o2e 
and the e2o configuration. (2) An overall odd number of 
lattice sites M (e.g. the o2e and the e2o configurations) 
shifts the filling factor in the leads away from 0.5 due to 
their finite size. A similar effect results from applying 
a gate voltage Vg 7^ 0, which imposes a problem to the 
extraction of the linear conductance. A possible solution 
is discussed in Sec. IV Bl 



Even-odd effect 



In [35| , a detailed analysis of finite size effects resulting 
from an even or odd number of lattice sites in the leads 
for a single-dot and for a three-dot structure with on-site 
interaction including the spin degree of freedom has been 
carried out. The behavior of the time dependence of the 
current resulting from the type of the lead (even or odd 
number of sites) has been traced back to the different 
magnetic moment of the system which is S^^^^^ — 1/2 for 
an overall odd number M of lattice sites and S^^^^y = for 
M being even. The reduction of the current in a situation 
where the leads both consist of an even number of sites 
(ene) as compared to the other possible combinations 
(one, ono) has been explained by the accumulation of 
spin on the structure in the first case corresponding to 
the effect of applying an external magnetic field. 

We already find parity effects in the time dependence 
of noninteracting spinless fermions in a system with a 
single-dot or a two-dot structure. Figs. [TTl [131 In the 
following we will trace the parity effects back to the 
level structure in the leads. The single particle levels 
ek of an uncoupled, noninteracting lead with Mi sites 
{i — L,R) arc given by = — 2i cos[7rfc/(Mi -I- 1)], 
k — 1 , . . . , Mi . The energy of a particle residing on a 
decoupled single dot structure (tc = 0) is simply given 
by the gate voltage Cd = V^, which is at the Fermi edge 
for Vg = {). For a decoupled n-dot structure one gets 
td,j = —2ts cos[7rj/ (n -I- 1)] -I- V^, j = 1, . . . , n. For an 
equal number of sites on both leads (as for example ene 
or ono) there is a twofold degeneracy of the single parti- 
cle lead levels which does not exist if Ml = Mr ± 1. In 
the degenerate case, single particle eigenfunctions can be 
constructed with a fully delocalized particle density while 
for Ml = Mr±1, the density profile of the single particle 
wave functions shows an alternating confinement of the 
particle on cither the left or the right lead The same holds 
true for the energy levels of the structure: if degenerate 
with a lead level, the single particle wave function can 
be distributed over the whole lead while it is localized 
on the structure otherwise. Therefore, in the ele case, 
the single-dot level is not degenerate with the lead levels 
when Cd = 0. As a result, a single particle occupying the 
dot level generates a sharp peak in the density profile (as 
well as the spin profile). For the olo case on the other 
hand, both leads have one energy level in the middle of 
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FIG. 11: Current through the contact link of a structure 
with two dots {ts = t), coupled to leads with a finite number 
of sites M and tc = 0.5t (compare also Fig. [TJ, at nomi- 
nal half filling N/M — 0.5 obtained from exact numerical 
diagonalization for bias voltage Vsd ~ 0.05t. The horizontal 
dotted lines represent the analytical result 7lb obtained from 
the Landauer-Biittiker approach. The current is measured 
on the left link to the structure. The time axis is normal- 
ized to the transit time Tr — Mh/{2t). Here, the focus is on 
finite size effects in the low voltage regime. We distinguish 
three cases: the system size is very small in case (a) where 
M = + X with a; = (29 lattice sites on the left and right 
which is an odd number in both cases o2o), x = 1 (now 30 
sites on the left which is an even number e2o), x = 2 (e2e) 
and a; = 3 (o2e). Here, the single particle level spacing Ae 
is much longer than Vsd, while the period of the Josephson 
oscillations Tj = 27rft/|VsD| is much bigger than the tran- 
sit time Tr. Case (b) shows an intermediate situation with 
M = 252 + x lattice sites. Here, Ae « Vsd and Tj f» Tr. A 
situation where Ae < Vsd and Tj < Tr is realized in case (c) 
with M = 1200 + X. For the e2o and the o2e case one has 
to do a density shift correction of the result since the total 
number of particles A'' / M/2, cf. Sec. IVB] 



the band, which together with the dot level generates 
a threefold degeneracy. For finite coupling tc > 0, the 
degeneracy of the lead levels and of the levels of the struc- 
ture with the lead levels gets lifted. The single particle 
wave functions must be divided equally on both leads, 
when A/l — A/r, while the alternating confinement is 
preserved for Ml = Mr ± 1 . Concerning the energy level 
of the dot, the threefold degeneracy in the uncoupled olo 
case results in two levels with strong localization on the 
dot, one lifted above the Fermi edge and one pushed be- 
low, and a third level with vanishing particle density on 



level occupation 
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FIG. 12: Initial occupation of the single particle energy levels 
in the non-interacting RLM (tc = 0.4t) at half filling. The 
number of lattice sites is M — Ml + Mr -|- 1 with the number 
of lattice sites in the left (right) lead Ml (Mr), (a) Ml -|- 
1 = Mr = 30. The alternating occupation can be traced 
back to the alternating localization of the single particle wave 
functions in either the left or the right lead, (b) A/l = Mr. = 
30. In the uncoupled case (tc ~ 0), the energy levels of the 
leads are degenerate. Therefore the energy levels can not be 
associated with only one lead. 



the dot, remaining on the Fermi edge. 

In a system with an odd number of lattice sites M and 
spinless electrons, half filling can not be realized strictly 
since N = AI/2 is not an integer. Adding spin shifts the 
particle number at half filling to A^ = M but leaves a 
total spin S^^^ — ±1/2, which will occupy the highest 
single particle level. Since for the doubly occupied levels 
the spin adds up to 0, the level at the Fermi edge de- 
termines the spin density profile which then explains the 
density peak on the dot in the ele case and the absence of 
a peak in the olo case. The time dependent behavior of 
the current can now be traced back to the single particle 
energy levels being confined in a single lead (fully dclo- 
calized) in the case of different numbers of lattice sites 
A/l = Mr ± 1 (equal number of lattice sites Ml = Mr). 
For the eno and one configurations, applying a bias volt- 
age as in Eq. ([7]) leads to an alternating occupation of 
the energy levels corresponding to the alternating con- 
finement of the single particle wave functions in the left 
or the right lead. In contrast we find an occupation num- 
ber of 1/2 in the energy range — Vsd/2 . . . Vsd/2 when 
Ml = Mr, corresponding to the fully delocalized single 
particle wave functions. We demonstrate this behavior 
for the non-interacting resonant level model (RLM) in 
Fig.IIl 

So far, we have a connection of the degeneracy of the 
single particle energy levels for the situation where the 
impurity is decoupled from the leads with the respec- 
tive class of the system (eno / one, ono, ene). The sit- 
uation changes when adding a constant local potential 
AV = AViiVL + AVr A^R to both, the initial and the 
time evolution Hamiltonian. To obtain the data of the 
dotted lines in Fig. [T3] we calculated the single particle 
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FIG. 13; Current through a structure coupled to two leads 
(mean value of left and right contact link) with an overall 
finite system size M at half filling obtained from exact di- 
agonalization. The figure demonstrates the influence of the 
number of lattice sites in the leads (even or odd) on the cur- 
rent for a bias voltage Vsd smaller than the single particle 
level spacing. The dotted lines represent a situation where 
an additional constant voltage is applied to both leads 

(a) or to the left lead (b), respectively. l\V 7^ results in a 
shift of the single particle levels in the uncoupled leads which 
can be used to "mimic" the different combinations of leads 
with an even or odd number of lattice sites, (a) M = 60 -I- a;, 
a; = (o2o), 1 (e2o), 2 (e2e) and 3 (o2e) where the num- 
ber of electrons is A*' = 30 for M = 60, 61 and A'^ = 31 for 
M = 62, 63. The dotted lines all together are generated us- 
ing a system with AI — 60 lattice sites, with AV 7^ 0. The 
different situations e2o and o2e can be recovered by changing 
the particle number from N = 30 to N = 31, cf. Sec. IV Bl 

(b) M = 61 + X, X = (ele), 1 (ole) and 2 (olo) where the 
particle number is fixed to A'' = 31. Here, the green (red) 
dotted line is generated from the ele (olo) system. 



energy levels for a system with an even (odd) number of 
lattice sites in the leads and then applied a relative shift 
of the lead levels with AVl = -AVr, G {e/4, e/2} for the 
two-dot structure and AVl G {ie/2}, AVr — for the 
single dot structure, where e is the energy gap to the first 
unoccupied energy level. This allows to change the level 
structure of a certain lead configuration in a way that it 
resembles one of the other configurations in the vicinity 
of the Fermi edge without changing the number of lattice 
sites in the leads. In Fig. [13] we see that the time depen- 
dent behavior of the system on the time scale T < Tr is 
only given by the structure of the single particle energy 
levels that contribute to the current, and the bias voltage 
Vsd, at least as long as we do not include interaction. We 
therefore conclude that ono as well as cnc configurations 
also can be used to study the I-V-characteristics in the 
low voltage regime. This may be interesting when inves- 
tigating structures with an even number of lattice sites 
on the structure, when the constraint N = AI/2 has to 
be fulfilled strictly. 
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FIG. 14: Current through a single impurity with an applied 
gate voltage — 0.21t for Vsd = 0.5t, coupled to two leads 
{tc = 0.3t), as a function of the system size. The analytic re- 
sult is obtained using the Landauer-Biittiker formula. While 
for different fiUings {N = M/2 and N = M/2 - 1) there is a 
systematic deviation from the analytic result, the interpola- 
tion results in a substantial improvement. The linear envelope 
is plotted to highlight the l/A'/-dependency of the finite size 
effects. For an explanation of the sinusoidal oscillations see 
also Fig. 1151 and the text. 



B. Density shift in the leads resulting from finite 
system size 

For the single resonant level model (RLM) the condi- 
tion of half filling is easily fulfilled by setting the particle 
number N = M/2 as long as the dot level resides in the 
middle of the band. Then the overall particle number 
density is n = 1/2 in the equilibrium case. This can 
change for diflFerent reasons: for example, for a model 
with two lattice sites in the structure and an overall odd 
number of lattice sites as discussed before half filling is 
not realisable, since AI/2 is not an integer. But even 
for the RLM, applying a gate voltage Vg ^ changes 
the particle number on the structure by AiVDot while 
changing the particle number per site in the leads by 
— AA^Dot/ {M — 1) which shifts the lead filling away from 
1/2 as long as the system size M is finite. In this section 
we will concentrate on the latter case. 

The impact on the current can be quite large, compare 
Figs. [HI [151 The total number of particles must therefore 
be corrected in such a way that A^Loads/(-^^ ~ 1) = 1/2 
where A^Leads = N — TVoot is the particle number in the 
leads. Thus an initial state has to be a mixture of 
states with different particle numbers I^'at) and |\E'Ar+i), 
or respectively, depending on the sign of AA^Dot 
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FIG. 15: Current through a single impurity with an apphed 
gate voltage Vg = 0.21t, coupled to two leads {tc = 0.3t), as 
a function of the voltage T^sd. The analytic result is obtained 
using the Landauer-Biittiker formula. The vertical lines rep- 
resent the single particle energies of a system with uncoupled 
leads {tc ~ 0.0); we find that the interpolated value of the 
current fits best with the analytical result if the bias voltage 
is chosen as the mean value of two neighboring energy levels 

(a) . However, this condition restricts the bias voltage to only 
a few values. The restriction can be circumvented by either 
increasing the number of lattice sites M or by using damped 
boundary conditions. The latter was used to obtain the values 

(b) without changing M - see section IVlI CI for discussion. 



For particle number conserving operators O the expecta- 
tion value reads 

(*i|d|*i) = \a\^{^N\d\^N) + |/3|'(«'jv±i|0|«'jv±i) 

(10) 

which leads to the condition 

|ap(*Jv|^Loads|*Jv) + 

+ l/5P(*JV±l|iVLeads|«'jV±l) = (H) 
W\' + m' = 1. (12) 

Since the current operator Ij also is particle number con- 
serving, the resulting time dependent current expectation 
value is an interpolation of the resuUs for N and for A^±l 
particles in the system 

Ij{T) = \a\^Ij{T;N) + {1 - \a\^)Ij{T;N ±1). (13) 

In Fig. [14] we show the dependency of the current 
through a single impurity coupled to two leads to the 
system size for different fillings N — M/2 as well as 
N = M/2 — 1, for a constant value of the bias voltage 
VsD and the gate voltage Vg. Furthermore we include the 
interpolated value, following the procedure described be- 
fore. We find that the interpolated results are centered 
around the analytic value, in contrast to the case with 
fixed particle number. However a distribution with an 
amplitude oc 1/M remains. A potential relation of the 
sinusoidal oscillations in the original data to the relative 



position of Vsd/2 to the single particle energy levels is 
illustrated in Fig. [151 Here, we show the current as a 
function of Vsd with Vg ^ 0, where we also apply the in- 
terpolation procedure. We compare the analytical result 
obtained using the Landauer-Biittiker approach with nu- 
merical data for the current through a single impurity 
coupled to two leads with a system size of M = 62 lat- 
tice sites in total. In order to interpolate the current as 
described before, Eq. (fT3)) . we simulated the time evolu- 
tion of the current expectation value with TV = 30 and 

= 31 particles in the system. In comparison to Fig. [Ml 
we conclude that one has to choose the system size in 
relation to the bias voltage carefully to get the desired 
relation of Vsd and the single particle levels. More pre- 
cisely, the data points (a), that fit nicely with the analytic 
curve, correspond to the interpolated current obtained 
for a bias voltage where Vsd/2 has been chosen as the 
mean value of two neighboring energy levels of the un- 
coupled (tc = 0) system. Another possibility is the use 
of damped boundary conditions to shift the single par- 
ticle levels, which yields the data points (b). This idea 
will be discussed in Section IVII CI 

A generalisation of this concept to systems with struc- 
tures of MDot > 1 sites with a corresponding number 
of energy levels is straightforward. A varying gate volt- 
age will change the occupation of the structure in a range 
A^Dot £ [0, AfDot] with a corresponding change of the par- 
ticle number in the leads. To get reliable results for the 
current at half filling in the leads it is then necessary 
to perform an interpolation of currents with appropriate 
particle numbers. Results for the linear conductance of 
a 7-site structure are discussed in the next section. 



VI. RESULTS FOR THE CONDUCTANCE 

Our result for the conductance through a single im- 
purity in Fig. [TBI is in excellent quantitative agreement 
with exact diagonalization results already for moderate 




FIG. 16: Current and differential conductance as function of 
applied potential through a single impurity with Vg = and 
half filled leads: N/M = 0.5. Circles (squares) show results 
for tc = 0.5t (0.35t). System size was M = 48 (M = 96) and 
Acut = 200 (400) states were kept in the DMRG. Lines are 
exact diagonalization results for M — 512. 
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FIG. 17: Differential conductance as a function of bias 
voltage through a 7 site nanostructure with nearest neigh- 
bor interaction. Parameters are tc ~ 0.5t, ts = 0.8t, and 
N/M=0.5. Squares (circles) denote weak (strong) interac- 
tion with U/ts = 1 (3) (here: Uc = 0.0). Lines are 
fits to a Lorentzian with an energy dependent self energy 
E = irjo + irjifj,^. Dashed lines: 771 — 0. System size is 
M = 144 {M = 192) and 600 (800) states were kept in the 
DMRG. 



system sizes and DMRG cutoffs. Accurate calculations 
for extended systems with interactions are more difficult, 
mainly for two reasons: 1.) The numerical effort required 
for our approach depends crucially on the time to reach a 
quasi-stationary state. For the single impurity, the quasi- 
stationary state is reached on a timescale proportional to 
the inverse of the width of the conductance resonance, 
4ifi/i^, in agreement with the result in Ref. [i^. In 
general, extended structures with interactions will take 
longer to reach a quasi-stationary state, and the time 
evolution has to be carried out to correspondingly longer 
times. 2.) In the adaptive td-DMRG, the truncation er- 
ror grows exponentially due to the continued application 
of the wave function projection, and causes the sudden 
onset of an exponentially growing error in the calculated 
time evolution after some time. This 'runaway' time is 
strongly dependent on the DMRG cutoff, and was first 
observed in an adaptive td-DMRG study of spin trans- 
port by Gobert ct al.[53|. To avoid these problems we 
resort to the full td-DMRG [s^l, which does not suffer 
from the runaway error. 

In Fig. [17] we show results for the first differential 
conductance peak of an interacting 7-site nanostructure. 
Gareful analysis of the data shows, that in order to re- 
produce the line shape accurately, one has to introduce 
an energy dependent self energy for U/ts = 3. Since 
the effect is small, we approximate it by a correction 
quadratic in the bias voltage difference fi = Vsb ~ Vpcuk- 
It is important to note that for the strongly interacting 
nanostructure, U/ts = 3, the conductance peaks are very 
well separated. Therefore the line shape does not over- 
lap with the neighboring peaks, and the fit is very robust. 
Performing the same analysis for a non-interacting nanos- 
tructure with a comparable resonance width, we obtain 
negligible corrections to ?/i in the self energy, indicating 
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FIG. 18: Transport through a non-interacting 7-site nanos- 
tructure with tc = 0-5t and ts ~ 0.8t. The energy levels of 
the nanostructure are indicated by dashed vertical lines, (a) 
Linear conductance for different A'^. The result after applying 
finite size corrections is shown as straight line (see text for 
details), (b) Number of fermions on the 7-site nanostructure. 
(c) Density p = {N - Noot)/{M - Muot) in the leads. System 
size is M = 96 and the number of states kept in the DMRG 
is iVcut = 400. 



that the change of the line shape is due to correlation 
effects. 

The linear conductance as a function of applied gate 
potential can be calculated in the same manner, if a 
sufficiently small external potential is used. We study 
the same 7-site nanostructure as before, with interaction 
U = 0, and use a bias voltage of Vsd = 2 • 10~^. For 
half filled leads, the result for the linear conductance cal- 
culated with a fixed number of fermions, N/M = 0.5, 
is qualitatively correct, but the conductance peaks are 
shifted to higher energies relative to the expected peak 
positions at the energy levels of the non-interacting sys- 
tem (Fig.[TH|). Varying the gate potential Vg increases the 
charge on the nanostructure by unity whenever an energy 
level of the nanostructure moves through the Fermi level 
[Fig. [THl (b)]. The density in the leads varies accordingly 
[Fig. [18] (c)] . Since the number of fermions in the system 
is restricted to integer values, direct calculation of the 
linear conductance at constant p is not possible and one 
must resort to interpolation. Using linear interpolation 
in p{N, Vg) for N = 44 ... 48 yields our final result for 
the linear conductance at half filling [Fig. [18] (a)]. The 
agreement in the peak positions is well within the ex- 
pected accuracy for a 96 site calculation. Our results for 
the conductance through an interacting extended nanos- 
tructure arc presented in Fig. 1191 The calculation for 
the weakly interacting system requires roughly the same 
numerical effort as the non-interacting system. In the 
strongly interacting case, where the nanostructure is now 
in the charge density wave regime, the time to reach a 
quasi-stationary state is longer, and a correspondingly 
larger system size was used in the calculation. In both 
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FIG. 19: Linear conductance through an interacting 7 site 
system with tc ~ O.bt and ts = 0.8t for weak (squares) and 
strong (circles) interaction. System size is M = 96 (M = 192) 
and 400 (600) states were kept in the DMRG. Finite size 
corrections have been included. Lines are guides to the eye. 
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cases we obtain peak heights for the central and first con- 
ductance resonance to within 1% of the conductance for 
a single channel. 



VII. EXPONENTIAL DAMPING 

In this section we want to study the effect and possi- 
ble applications of damped boundary conditions (DBC). 
DBC have been introduced [H, [s^l in order to reduce 
finite size effects. Here we would like to reduce the lim- 
itations rising from the finite transit time Tr and the 
Josephson wiggling which especially in the low voltage 
regime and with an applied gate voltage spoils the accu- 
racy of current measurements. We have already seen how 
to profit from the voltage dependency of the finite size 
wiggling by using a fit procedure which allows for the cal- 
culation of current -voltage characteristics even with an 
applied gate voltage. We now want to discuss the pos- 
sibility of combining the fit procedure with DBC, where 
the damping effectively increases the system size. Fur- 
thermore we want to use DBC to adjust the single parti- 
cle energy levels in order to increase the resolution with 
respect to Vsd when Vg ^ 0, cf. Fig. [T5l 



A. Estimate for Transit Time in a system with half 

filling 

In Fig.[5D]we show the time dependent current through 
a single impurity with Vg = 0, including the initial tran- 
sient regime as well as the finite size reflections for dif- 
ferent values of the bias voltage Vsu- We compare two 
different system sizes with M = 120 and M = 240 
lattice sites, and also apply exponentially DBC in or- 
der to demonstrate the effectively increased system size. 
The hopping matrix element is damped towards the 
boundaries of the system using a damping constant A 
as sketched in Fig. [21 over a range of Ma lattice sites. 
The total number of lattice sites is left unchanged (here: 
M = 120). We find an enhanced size of the current 
plateaus, however, the damping can also lead to an early 



FIG. 20: Time dependent current through a single impu- 
rity with tc = 0-St at nominal half filling N/M = 0.5 ob- 
tained from exact numerical diagonalization for different bias 
voltages Vsd and different damping conditions. For small 
bias voltage, finite size reflections from hard wall boundary 
conditions (HWBC, a) can be suppressed significantly using 
damped boundary conditions (DBC). Using an exponential 
damping with A^^^^ = 0.93, M = 120 and A/a = 50 (b) 
yields a plateau of constant current for Vsd = 0.4t consider- 
ably bigger than in the undamped case. However, the current 
plateau starts dropping before the estimated transit time ac- 
cording to Eq. (|15|) is reached (here: Tr « 670), which gets 
even more pronounced when increasing the bias voltage. Re- 
ducing the damping (c, d) can lead to good agreement with 
the estimate (rR(c) « 178, Tn{d) « 123). 



breakdown of the current. 

As an estimate for the transit time of a wave packet 
traveling in undamped leads of size M one can use the 
Fermi velocity vp = 2t/h which leads to 



M 

Vf 



Mh 
'2F' 



(14) 



Assuming a local Fermi velocity vp{x) = 2t{x)/h in 
damped leads with damping A > 1 leads to an expression 
of the form 



Mh 

~2r 



1 - 



2Ma\ 



2h 



M 



') HnA 



where Ma is the size of the damped leads. Eq. (|T5|l can 
then be used to estimate an effective system size 



M, 



off 



M - 2Ma + r—r lA'"'^/'^ - 1) 
m A ^ ' 



(16) 



This is in agreement with the results for the pseudo- 
steady current found for the noninteracting case, Fig.[20l 
For a more quantitative check of the formula we compare 
the transit time, extracted from a current measurement, 
to the estimate given by Eq. ((T5)) [Fig.[2T]. We therefore 
use two different criteria: (a) the time t'^'^ where /(T) 
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FIG. 21: Test for the transit time estimate Tr of the current 
through a single impurity at half filling, Eqns. (|14lll5p . where 
the black dotted line is the undamped case. All values are 
plotted as functions of the damped lead size Ma. The small 
plots at the top show the single particle level density for the 
energy given by the bias voltage, in units of the level density 
for the undamped case. (See text for details) 



becomes negative at the end of the first plateau (crosses) , 

and (b) the time where the current changes sign af- 
ter one round trip (squares). The black dotted lines show 
and for the undamped case. For values of A~^/^ 
close to 1 we find that the estimate is well fulfilled over 



a wide range of values of A/a for both (a) and (b) even 
for big bias voltage. The slight growth of Tp'^'^^'/TR is 
assumed to be caused by the different Fermi velocity of 
excitations for |Vsd| > 0. However, the estimate tends 
to be totally wrong even for small bias voltage and small 
values of Ma if A~^/^ becomes too small. The small 
plots at the top show the relative single particle level 
density. As expected, cf. Fig. [22l the level density grows 
with Ma until a maximum is reached where the position 
of the maximum is determined by the bias voltage. It 
can clearly be seen that the position of the maximum in 

combination with the values of /T^, gives a strong in- 
dication if a current plateau is still well defined for a time 
scale given by the estimate of Tr, since t!^^ /Tyi ~ 1 for 
values of Ma on the left side of the maximum of the single 
particle level density. In comparison, (b) is a weak crite- 
rion since for strong damping the current plateau starts 
decaying for times much smaller than Tr, cf. Fig. [201 
In Fig. [2H we show the single particle energy levels of a 
system with M = 120 lattice sites with a single impu- 
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FIG. 22: Level discretisation in a finite system (M = 120) 
with a single impurity, coupled to leads (tc = 0.3) as function 
of the damping rate A^^'''^ (a, b), as well as function of the 
size Ma of the damped leads (c). The damping lead size is set 
to (a) Ma = 30 and (b) Ma = 50, while for (c) the damping 
rate is set to A~^'^^ = 0.98. The implementation of damped 
leads in combination with leads described by a uniform tight 
binding chain can be used to increase the level density in the 
vicinity of the fermi edge while allowing for direct access to 
real space quantities like the current at a specific lattice site, 
as e.g., the impurity. However, this approach is only useful for 
the calculation of current in a limited voltage window, since 
in the high voltage regime also energy levels at the band edge 
get occupied, where the level spacing is significantly increased 
with A and Ma. 



rity, as function of the damping constant A^^/^ as well 
as function of the size of the damped leads Ma- The plot 
demonstrates the growth of the level density on the scale 
-^jiich in conjunction with Fig. [21] allows for an 
estimate of the maximum value of Vsd up to which a 
current plateau can be expected in a system with DEC. 



B. Fit Procedure 

As already mentioned in Sec. [Vl the fitting procedure 
gets unreliable when the oscillation time Tj substantially 
exceeds the time range Ts . . . Tr. We therefore now want 
to demonstrate how to use the estimate for the tran- 
sit time in order to implement damping conditions to 
sufficiently increase the effective system size, enforcing 
Tj ~ Tr — Tg. As an example, we simulate the time 
evolution of a system with M lattice sites and a single, 
non-interacting impurity with Vg = 0, and apply a small 
bias voltage Vgo > 0. An effective transit time k, Tj 
can be obtained using DEC, according to Eqns. ([T51 [TB)) . 
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FIG. 23: Current through a single impurity with tc ~ 0-3t 
and Vg = 0. The time axis is normalized to the oscillation 
period Tj = 2'n-h/Vso, with (a) Vsd = 0.02t and (b) Vsd = 
O.lt. The analytic results are computed using the Landauer- 
Biittiker formula. For Vsd = 0.02t (a), the oscillation period 
is Tj = 314:h/t. To obtain a current plateau containing at 
least one Josephson oscillation one has to simulate the time 
evolution of a system with AI > 630, which is very hard on 
present days computers when interaction is included. Here, 
we apply DBC on a system with M = 96 (M = 192) to 
effectively increase the system size using (i) A ~ 0.903, Ma = 
32 (ii, A 0.969, Ma = 84). Accidentally, the fit value agrees 
with the analytic value nearly perfectly for configuration (i). 
For Vsd = O.lt (b), Tj = 63Vi ^ M > 126. The damping 
conditions are characterized by (iii) A « 0.93, Ma ~ 20 and 
(iv) A « 0.900, Ma ~ 20, respectively. In both cases one is 
able to extract a current via the fit procedure although M <C 
2tTj/h. However, the most reliable results can be obtained 
by inceasing M, c.f. M = 180 in (b). 

The result is presented in Fig. I23[ where we show 
the time dependent current through one of the contaet 
hnks of a single impurity for different damping condi- 
tions and two different values of Vsd- Again, we fit 
I + Ij cos(VsDr + 1^) to the oscillating part of current 
expectation value. The extracted current / for the calcu- 
lations including DBC fits with the analytic result with 
an accuracy of ^ 1% which is of the same order of mag- 
nitude as compared to the mean value of the very small 
plateau regime that can be found for the system with 
HWBC. This leads us to the conclusion, that DBC can be 
used to obtain a first guess while for high precision mea- 
surements, HWBC with an increased system size have to 
be implemented. 



C. Correction of the single particle energy levels 
using DBC 

In Section |VB] we found that the effects resulting from 
a finite density shift in the leads when applying a gate 
voltage can be significantly suppressed when extracting 
the current only for certain values of Vsd determined by 
the single particle level spacing. Since these finite size 
effects particularly arise in the middle of the band where 
the density of single particle levels is the lowest - and 
where the current has to be extracted for the calculation 
of the linear conductance - one would like to shift the 
levels towards the center of the band somehow. This 
can be achieved by increasing the number of lattice sites 
which also increases the numerical effort. 

Applying DBC also results in a shift of the single par- 
ticle energy levels in the leads towards the center of the 
band, cf. Fig. We therefore state the question if the 
criterion formulated in Sec. IV Bl still holds true for DBC. 
The result is shown in Fig. [TH To obtain the additional 
data points (b) we used damping conditions with values 
of A-i/2 = 0.91 . . . 0.98 and Ma = 15, 20, 23. We calcu- 
lated the single particle energy levels for the decoupled 
leads and then obtained the current for values of the bias 
voltage with Vsd/2 in the middle of two neighboring en- 
ergy levels. To increase the resolution for the high voltage 
regime only moderate damping conditions are required 
(A-i/2 = 0.98, Ma = 15,20), while strong damping is 
imposed to get high resolution in the low voltage regime. 
For Vsd approaching the band edge, however, DBC have 
to be avoided for the reasons discussed above. 



VIII. CONCLUSIONS 

We have reviewed the concept of extracting the finite 
bias and linear conductance from real time evolution cal- 
culations in finite systems. Very accurate quantitative 
results are possible, as long as finite size effects are taken 
into account. Our results for the linear conductance com- 
pare favorably both in accuracy and computational effort 
with the DMRG evaluation of the Kubo formula [l3| . 
Calculations of strongly interacting systems show corre- 
lation induced corrections to the resonance line shape. 
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